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ABSTRACT 


A  method  for  accurate  numerical  integration  is  proposed  and 
illustrated  by  a  variety  of  examples.  The  method  depends  on  an 
ability  to  efficiently  evaluate  the  weights  and  abscissas  of 
Gaussian  quadrature  formulas,  and  software  to  achieve  this  pur¬ 
pose  is  also  presented. 
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INTRODUCTION 


A  considerable  number  of  problems  of  physical  and  engineer¬ 
ing  interest  require  as  part  of  their  solution  the  accurate 
evaluation  of  definite  integrals  of  the  form 


f (x) dx. 


[*] 


where  R  is  a  closed  region  of  one  or  more  dimensions.  Typical 
cases  have  R  =  a  compact  or  infinite  interval  in  one  dimension, 
and  R  =  a  rectangular  or  spherical  solid  in  higher  dimensions. 

In  applications,  the  integral  [*]  may  yield  percentage  points 
of  a  probability  distribution ,  the  autocorrelation  of  a  signal, 
a  one  or  two  dimensional  convolution  of  a  signal  with  a  filter 
response  function,  the  intensity  of  graybody  radiation  emitted 
into  a  particular  waveband  (Planck  integral) ,  certain  thermo¬ 
dynamic  lattice  sums,  etc. 

Usually  it  is  desired  that  the  computed  value  of  [*]  be 
accurate  to  so  many  digits  (prescribed  relative  error)  or  to 
so  many  decimal  digits  (prescribed  absolute  error) .  The  only 
way  that  a  user  can  be  certain  that  this  requirement  is  satisfied 
is  to  have  available  an  error  formula  for  the  integration  rule 
being  employed,  and  to  verify  that  the  integrand  and  the  order 
of  the  rule  are  such  that  the  error  term  is  sufficiently  small. 
(In  making  this  last  statement  we  are  assuming  negligible  round- 
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off  error.)  Now  typically  the  classical  error  term  (  for  inte¬ 
gration  over  intervals)  involves  a  small  coefficient  multiplied 
by  a  high  order  derivative  evaluated  at  an  unknown  point  in  the 
interval  of  integration.  Hence  in  practice  the  maximum  modulus 
of  this  derivative  must  be  computed  or  at  least  bounded.  And 
in  practice  this  task  is  usually  impossible,  due  to  the  complexity 
of  the  intergrand. 

One  way  out  of  this  difficulty  (for  the  case  where  R  is  a 
compact  interval)  is  available  when  the  integrand  f  extends  to 
an  analytic  function  on  some  neighborhood  of  R.  Complex  variable 
methods  (related  to  the  Cauchy  integral  formula)  then  often  per¬ 
mit  bounds  on  the  error  in  terms  only  of  the  values  of  f,  and  not 
of  its  derivatives.  Furthermore,  the  error  term  will  tend  ex¬ 
ponentially  to  zero,  as  the  order  of  the  integration  rule  in¬ 
creases.  Unfortunately,  it  is  not  always  the  case  that  our  in¬ 
tegrand  is  analytic. 

This  note  will  present  a  procedure  for  resolving  the  problem 
of  fast  accurate  quadrature  for  integrands  which  are  reasonably 
smooth  but  not  necessarily  analytic.  In  the  absence  of  smooth¬ 
ness  (i.e.,  if  there  are  discontinuities  in  the  first  derivative, 
or  in  one  of  the  first  partial  derivatives,  if  the  integrand  is 
a  function  of  several  variables) ,  our  method  will  not  have  favor¬ 
able  convergence  properties,  and  should  not  be  employed.  In  the 
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case  of  integration  over  a  compact  interval,  one  would  then  be 
better  advised  to  utilize  a  trapezoidal-type  rule,  such  as 
(adaptive)  Romberg  integration  [1,  §  6.3].  Software  for  such  a 
procedure  is  available  in  the  International  Mathematical  and 
Statistical  Library  (IMSL)  under  the  title  DCADRE . 

In  brief,  our  method  is  based  on  observation  of  the  con¬ 
vergence  of  a  sequence  of  Gauss-type  quadratures  with  appropriate 
weight  function.  For  this  to  be  practical  it  is  necessary  to 
have  available  a  fast  procedure  for  generating  the  abscissas  and 
coefficients  of  a  Gaussian  formula  of  a  given  order.  Such  a 
procudure  is  given  below;  it  is  based  on  suggestions  of  Wilf 
[2]  and  Golub-Welsch  [3]  to  convert  this  problem  to  an  eigen- 
analysis  of  a  certain  (symmetric)  matrix. 

The  original  motivation  for  this  work  came  from  the  (two- 
dimensional)  problem  of  accurately  computing  the  overall  output 
of  an  optical  detector  [4,  §  6. 2. 2. 1.1].  This  output  is  the  con¬ 
volution  of  the  detector  response  function  and  the  irradiance 
function  H.  In  the  important  special  case  of  a  rectangular 
detector  this  convolution  reduces  to  the  integration  of  H  over 
the  (rectangular)  detector  surface.  Our  method,  described  below, 
was  proposed  as  an  alternative  to  the  use  of  fast  Fourier  trans¬ 
form  techniques,  with  their  arbitrary  sampling  decisions  and 
apparent  lack  of  suitable  error  theory. 
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II.  REVIEW  OF  GMJSSIAN  QUADRATURE  THEORY 


This  section  contains  a  brief  exposition  of  convergence  and 
error  analysis  for  Gaussian  integration  formulas.  No  attempt 
has  been  made  to  write  a  small  text  on  the  subject!  For  further 
background  the  sources  [1,  5,  6,  7]  may  be  consulted.  In  general 
the  recent  monograph  [1]  is  an  excellent  source  of  information 
about  all  aspects  of  the  quadrature  problem. 


2.1  Basic  Theor\ 


Let  I  be  an  interval  of  real  numbers  and  let  w  be  a  weight 

function  on  I,  that  is,  oj(x)>0  and  y*xnw(x)dx  is  finite  for 

I 

n  =  0,1,2,**«.  Then  there  exists  a  sequence  {pn (x) :  n  =  0,1,2, 
. . . }  of  polynomials  of  degree  n  with  are  orthonormal  on  I  with 
respect  to  w,  in  the  sense  that 


ra 


Pn(x) 


a)  (x)  dx  = 


(2.1) 


for  all  m,n,  =  0,1,2,  ••>•.  These  pn  are  unique  up  to  a  sign,  and 
we  shall  adjust  this  sign  so  that  the  coefficient  kn  of  xn  in 
Pn  is  positive.  The  pn  can  be  computed  either  by  applying  the 
Gram-Schmidt  procedure  to  the  monomials  {xn}  or  by  use  of  a 
recurrence  relation 


(2.2) 


pn(x)  =  (anx  +  bn)  pn-1(x)  -  cn  Pn_2(x), 


where  an»cn  +  Of  p_3(x)  =  °/  Pq  (x)  =  (x)  dx)  .  This  re¬ 

currence  will  be  very  important  to  us  below.  Let  us  now  just 

note  one  special  relation  between  the  coefficients  a  and  c 

n  n 

that  follows  from  the  normality  of  the  pn»  Namely,  by  multiply¬ 
ing  both  sides  of  (2.2)  by  pnw  and  by  Pn_2w  an<3  integrating  over 
I,  we  obtain 


1 


0 


y*xpn_1(x)  pn_2  (x)  to  (x)  dx. 


i  /  xp  . 
n  J  *n-l 


a_  /  xp_  ,  (x)  pn_2 (x)u  (x)dx  -  cn, 
I 


whence 


c  =  a  /a  . 
n  n  n-1 


(2.3) 


Now  consider  the  n  x  n  tridiagonal  symmetric  (Jacobi)  matrix 


A  = 


--Vai 

c2/a2 

0 

0 

0 

C2/a2 

~^2^a2 

c3/a3 

0 

0 

0 

c3/a3 

-b3/a3 

c4/a4  ** 

0 

0 

0 

•  •  • 

c  /a 
n  n 

'  Van  - 

Using  (2.3)  the  reccurrence  relation  (2.2)  can  be  rewritten  in 
the  matrix-vector  format 


AP  (x) 


xP  ( x ) 


+ 


Cn+1 

an+l 


(2.4) 


T 

where  P(x)  =  (pQ(x),  p-^x),...,  p  ^(x))  .  From  (2.4)  we  con¬ 
clude  that  a  real  number  xq  is  a  root  of  p^  exactly  when  it  is 
an  eigenvalue  of  A.  [2,3], 

The  reason  for  our  interest  in  the  roots  of  p  is  of  course 

rn 

that  they  comprise  the  abscissas  of  the  corresponding  Gauss 

fch 

quadrature  rules.  The  n  Gauss  quadrature  formula  is  a  (pos¬ 
itive)  linear  functional  whose  value  at  a  continuous  function 
f  is 


where 


Gn(f;u)) 


(2.5) 


0  <  w. 


(n)  _  _  n+1 


n 


(P  (x  ^  ) 
'  n+llxk  ] 


n 


(x. 


(n) 


-1 


)) 


l<k<n. 


With  this  choice  of  weights  { (n )  }  an(j  abscissas  {x^n^  }  it  is 
well  known  [1]  (see  also  §  3.2  below)  that 
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J  f(x)w(x)dx  =  G  (f;w) 
~  n 


+  f (2n) (g) 

(2n) 1  k  2 
n 


a<  r<b 


(2.6) 


provided  that  a,b  are  finite  and  f  has  2n  continuous  derivatives 
on  [a,b].  A  particular  import  of  (2.6)  is  that  the  Gauss  rule 
of  order  n  exactly  integrates  every  polynomial  of  degree  <  2n  -  1, 
This  is  the  famous  optimality  property  of  Gauss  rules:  they  are 
exact  for  polynomials  of  as  large  a  degree  as  is  possible  with 
a  formula  of  the  type  (2.5)  (there  are  only  2n  free  parameters 
in  (2.5)). 


2.2  Error  Analysis 


Let  En(f)  denote  the  error  in  uniform  approximation  on 
[a,b]  to  the  continuous  function  f  by  polynomials  of  degree 
<_  n.  If  p  is  such  a  polynomial  then  the  error  in  ntl1  order 
Gauss  quadrature,  EGn(f;oj),  satisfies 


EGn(f?u))  | 


'/ 


b 

f  (x)  03  (x)  dx 


-  Gn(f;cu) 


< 


b 

I J  <f  (x)-p(x)  )u>(x)dx|  +  |Gn(f-p;w) 
a 


< 


max  | f (x)-p (x)  j 
a<x<b 
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This  being  true  for  all  polynomials  of  degree  <  2n-l,  we 
conclude  that 


|EGn(f;a))  |  <  2 \j 
a 


to  (x)  dx 


E2n-l(f) 


(2.8) 


The  advantage  of  (2.8)  over  (2.6)  is  first,  that  there  is  no 
derivative  assumption  made  about  f,  and  second,  that  upper  bounds 
for  the  error  can  always  be  obtained  by-  use  of  (2.7).  Of  course 
the  best  theoretical  estimates  of  ER(f)  depend  on  derivative  in¬ 
formation  about  f  ("Jackson's  theorems",  [6]). 

In  the  most  favorable  situation  for  Gaussian  quadrature  on 
a  compact  interval  [a,b],  the  integrand  f  is  (the  restriction  of) 
an  analytic  function  on  [a,bj.  Then  it  is  known  ("Bernstein's 
theorem",  [6])  that  the  quantities  En(f)  decrease  very  rapidly; 
precisely. 


En  (f )  1  Kq11, 


for  positive  constants  K  and  q,  with  q  <  1.  Hence,  by  virtue 
of  (2.8)  the  Gauss  quadratures  of  f  converge  at  essentially  a 
geometric  rate  ("r-linear  convergence")  to  the  true  value  of 


f  b 

J  f(x)ui(x)dx.  (There  is  in  fact  an  extensive  literature  on  error 
a 

estimates  for  Gauss  quadratures  of  analytic  functions,  giving 
more  precise  versions  of  this  fact;  see  [8]  for  a  recent  con¬ 
tribution.  ) 

The  last  two  paragraphs  strongly  suggest  that  Gaussian 
quadrature  will  be  most  effective  for  integrands  which  exhibit 
polynomial  behavior,  and  this  is  certainly  true  in  a  practical 
sense,  when  computation  time  is  an  issue.  But,  in  fact,  the 
estimate  (2.8)  and  the  Weiers trass  approximation  theorem  show 
that 

lim  G 
n->°° 

holds  true  whenever  f  is  a  continous  function  on  [a,b].  Indeed, 
(2.9)  is  valid  whenever  fto  is  Riemann  integrable  on  [a,b]  [9], 

2.3  Double  Integrals 

Let  us  now  consider  the  problem  of  numerically  evaluating  a 
multiple  integral  over  a  rectangular  region  R  =  [c,d]  x  [a,b]. 

The  ensuing  discussion  applies  as  well  to  solid  rectangular 
regions  in  3  or  more  dimensions,  and  integrals  over  other  types 
of  regions  can  frequently  be  reduced  to  integrals  over  rectangular 
ones  by  means  of  appropriate  transformations  of  coordinates 
(1,  §  5.4]. 
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n 


(f  ;o>) 


-I 


f  (x)  a)  (x)  dx. 


(2.9) 


A  standard  approach  to  evaluation  of 


d 

f  (x,y)w  (x)lu  (y)dx  dy 


(2.10) 


is  to  apply  one-dimensional  quadrature  formulas  to  the  functions 
f(x,*)  and  f(*,y).  The  resulting  formula  is  known  as  a  product 
rule  [1,  §  5.6,  10,  §  2.3]  and  in  the  case  that  each  one¬ 
dimensional  formula  is  G  (•;u).  we  can  obtain  the  product  Gauss 

n 

formula 


n 


n 


G(2)  (f 
n 


,  V  V  (n)  (n)  (n) 

;w)  =?.  Zi  w  w.  f(x.  , 

i^l  3=1  1  3  1 


x<nl ) 


(2.11) 


as  an  approximation  to  (2.10).  Although  the  theory  of  multivariate 
polynomial  approximation  is  not  as  well  developed  as  the  single 
variable  case,  there  is  still  a  Weierstrass  theorem.  Hence, 
since  (2.11)  will  integrate  any  polynomial  p(x,y)  exactly  provided 
that  n  is  sufficiently  large,  we  may  conclude  that 


lim  G^n^  (f  ;w) 
n-*-°° 


d 

f  (x,y)u>  (x)  w  (y)dxdy. 


(2.12) 


for  all  continuous  functions  on  R. 

What  about  the  validity  of  (2.12)  when  the  integrand  in 
(2.10)  is  merely  Riemann  integrable?  The  difficulty  here  is  that 
this  assumption  need  not  imply  that  the  sections  f(x,y)co(x)  and 


k—.  * 


f(x,y)t»)(y)  be  integrable  on  [c,d]  and  [a,b],  respectively.  Hence, 
in  general  there  will  be  no  reason  to  expect  the  product  Gauss 
rule  to  converge.  But  suppose  that  the  integral  (2.10)  is  well 
enough  behaved  to  obey  the  Fubini  iterated  integrals  formula: 


f(x,y)u(x)oi(y)  dxdy 


a  c 


d 

f (x,y)w (x)dx)w (y) dy . 


(2.13) 


Denoting  the  inner  integral  in  (2.13)  by  f(y),  we  can  compute 

Gn(f;*)  for  a  sufficiently  large  n;  then  we  can  approximately 

/  „  \ 

compute  each  value  f (x^  ')  for  k  =  l,...,n  by  Gaussian  formulas 
/  n  \ 

G  (f  (ui,x1:  ))  for  sufficiently  large  m.  In  this  way  we  can 

III  K 

conclude  that 

n  m  b  d 

lim  lim  ^2  ^2  cj.oi.  f(xfm}  )  =f  f  f  (x,y)  uj  (x)  w  (y)  dxdy . 

n-*°°  m-><»  i=l  j=l  3  J  a  J  c 

This  convergence  result  is  clearly  less  satisfactory  than  (2.12) 
although  it  could  be  used  in  practice  provided  (2.13)  could  be 
verified. 


2.4  Infinite  Interval 

Finally  we  consider  the  matter  of  the  validity  of  the  con¬ 
vergence  formula  (2.9)  when  b  =  °°.  This  is  not  covered  by  the 
earlier  analysis  on  account  of  the  non-compactness  of  the  in- 
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terval  [a,“).  However,  it  is  shown  in  [5,  §  III.l]  that  for  a 
class  of  positive  quadrature  rules  which  includes  the  important 
special  case  of  the  Gauss-Laguerre  formula  with  weight  function 
to  (x; a)  =  xae  x  (a>-l)  ,  convergence  is  obtained  for  every  function 
f  such  the  f  ( • ) ci)  ( • ; a)  is  Riemann  integrable  on  [a,“) ,  and  f  has 
at  most  polynomial  growth  at  The  proof  depends  on  being  able 

to  "squeeze"  f  between  two  polynomials:  P1£f£P2»  such  that 
f  Tp2(x)  -  P-^  (x) ) a)  (x;a)dx  is  arbitrarily  small. 

III.  THE  PROPOSED  METHOD 

Armed  with  this  theoretical  underpinning  of  convergence 
theory  we  can  in  practice  proceed  to  evaluate  an  integral  of  the 
form  [*]  (Introduction)  by  simply  observing  the  behavior  of  a 
sequence  of  approximations  Gn(f;u)):  n=l,  2,...  ,  and  terminating 
the  process  when  |Gn+1(f;w)  -  G^  ( f  y a> )  |  <  e,  for  some  preassigned 
tolerance  e.  There  are  several  caveats  to  be  mentioned  here. 
First,  this  termination  criterion  cannot  conclusively  prove  that 
G^  ( f ; a) )  is  adequately  close  to  [*]  (cf.  the  discussion  in  [1,  § 
6.1]).  It  is  always  possible  to  construct  ad  hoc  counterexamples. 
Nevertheless,  this  does  not  seem  to  be  a  difficulty  in  practice, 
although  it  would  be  safer  to  require  that  the  termination 
criterion  above  hold  for  2  or  3  successive  values  of  n.  Second, 
there  is  the  matter  of  choice  of  weight  function,  especially  when 
one  is  not  obviously  present  to  begin  with.  Third,  there  is  the 


need  to  be  able  to  efficiently  generate  the  weights  and  abscissas 
necessary  to  fully  specify  the  value  Gn(f;w)  for  various  values  of 
n.  Finally,  this  method  should,  as  already  noted,  only  be  applied 
in  cases  where  the  integrand  f  is  reasonably  smooth;  otherwise 
convergence  is  likely  to  be  too  slow  to  be  useful. 

The  choice  of  weight  function  will  be  discussed  in  the 
following  sections.  As  to  the  third  point  just  raised,  the  key 
here  is  the  use  of  the  eigenvalue  problem  shown  in  formula  (2.4). 
By  computing  the  eigenvalues  of  the  Jacobi  matrix  A  we  obtain  the 
abscissas  of  the  corresponding  Gauss  formula.  Further  it  is 
known  that 

1  =  c^n)  P(x<<  P(x^n)), 
where  P(x)  was  defined  just  after  (2.4).  Hence 


is  a  normalized  eigenvector  of  A  corresponding  to  the  eigenvalue 
x^n) .  By  considering  just  the  first  component  of  P(x^n^)  we  can 
conclude  that 

=  J  a)  (x)dx  •  (Q,[n)  )1. 
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In  other  words,  the  various  weights  associated  with  the 

Gauss  formula  Gn  (•;<*))  can  be  obtained  as  the  constant  (x)  dx 
times  the  square  of  the  first  component  of  a  system  of  n 
orthonormal  eigenvectors  associated  with  the  matrix  A. 

FORTRAN  codes  illustrating  this  procedure  for  several  weight 
functions  and  numerical  illustrations  of  the  whole  method  will 
be  given  in  the  ensuing  sections,  along  with  further  discussion. 

It  is  important  to  appreciate  the  ease  and  accuracy  with  which 

the  weights  and  abscissas  can  be  computed  in  this  manner.  Earlier 
efforts  were  much  slower.  They  typically  involved  the  use  of 
some  root  finding  scheme,  such  as  Newton's  method,  applied 
directly  to  locate  each  zero  of  the  relevant  polynomial.  All  the 
usual  attendant  difficulties  arose  here,  such  as  the  need  for  a 
sufficiently  accurate  starting  value.  This  was  particularly  a 
problem  in  trying  to  locate  the  zeros  of  the  generalized  Gauss- 
Laguerre  polynomials  (corresponding  to  the  weight  function 
oj  (x;a)  =  xae_X,  a  >  -1,  on  [ 0 , °°) )  ,  as  these  tend  to  drift  off  to 
infinity  as  either  a  or  n  increase.  See,  for  example,  the  dis¬ 
cussion  in  [11]. 


IV.  GAUSS-CHEBYSHEV  QUADRATURE 


We  first  look  at  the  case  where  the  weight  function  u (x)  = 

2  _35 

(1-x  )  on  [-1,  1].  Thus  the  resulting  quadrature  formula  (2.5) 
will  approximate 


f  f (x) dx 
-1 


The  orthonormal  polynonials  associated  with  this  u>  are  the  well- 
known  Chebyshev  polynomials: 


Pn(x) 


cos (n  cos 


n>l 


P0(x)  = 


In  this  case  we  may  clearly  write  down  the  roots  of  p  (  =  the 
abscissas  of  the  corresponding  Gauss  formula)  by  inspection: 


x^n^  =  cos  (~2n~  1T)  r  l<k<n. 

Further,  although  less  obvious,  the  associated  weights  =  ir/n, 
l<k£n,  independently  of  k  [6,  v.2,  ch.5].  It  is  interesting  that 
this  is  the  only  case  where  this  property  of  the  weights  holds 
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*!■  a*#-  Hr 


true  (assuming  the  interval  is  fixed  as  [-1,  1  ];  this  is  Posse's 
theorm  [6,  v.2,  Ch.6]). 

Combining  all  this  information  the  general  formula  (2.6)  here 
appears  as 


f  f  (x)dx 
-i  yi-x2 


X] f  (c°s  < 

n  k=l 


2k-l 

2n 


it)  ) 


+ 


2tt  f  (2n)  (g) 
4n(2n)  i 


(3.2) 


for  some  |£|  <  1,  provided  f  e  C^n[-1,  1].  This  formula  is  re¬ 
markable  because  the  weights  and  abscissas  appear  in  closed  form; 
the  eigenanalysis  of  the  preceding  section  is  not  needed  here  for 
their  computation. 


4.1  Program  and  Examples 


Exhibit  1  gives  a  short  FORTRAN  program  which  enables  a  user  to 
calculate  an  ntl1  order  Gauss-Chebyshev  approximation  to  any  in¬ 
tegral  of  the  form 


(3.3) 


the  programs  accepts  a,  b,  and  n  as  inputs  and  calculates  the 
corresponding  quadrature  (after  an  appropriate  change  of  variable) 
Because  in  general  the  integral  (3.3)  does  not  contain  the  weight 
function  w  as  a  factor,  it  is  introduced  artifically  in  the  FUNC- 
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c 


1 


25 

2 


5 


ie 

12 

20 


PROGRAM  TO  DO  ONE  DIM.  GAUSS-CHEBVSHEU  QUADRATURE 


COMMON  A,  B 
URITE (6, 1 ) 

FORMAT! IX, 'ONE  DIN.  GAUSS-CHEBVSHEU  QUADRATURE 
MX, ' INPUT  LOWER  AND  UPPER  END  POINTS') 

RE AD <5, 1 )  A  , B 
URITEI6.2) 

FORMAT! IX. 'INPUT  ORDER  OF  G-CHB.  FORMULA') 
READ15,*)  N 

CALL  TIMES(Dl.Tl.IU.IT) 

PI2  -  AS1NI1.  ) 

ONF  ■  0. 

DO  5  K-l.N 

XICN  -  COSl  <2*K-1  )*PI2/N) 

QNF  •  QNF  ♦  F<XKN> 

CONTINUE 

QNF  •  QNF*<B-A5*PI2'-N 
CALL  TIMEStD2,T2,IU2,IT2> 

ETIM  »  FLOAT! 102- IU >*.000013 
URITECS.10)  N,  QNF 

F0RMATI/1X, 'G-CHB.  QUAD' F ; ' . 13. ' )  *  .F12.6) 

URITE<6, 12 )  ETIM 

F0RMAT(/1<, ‘EXECUTION  TINE  ■  '.F7.3) 

URITE ( 6,20 ) 

FORMAT' /IX, 'RESTART?  INPUT  0(NO»  OR  UVES>'> 

READ! 5, *  )  HP 

IF  (NP  .EQ.  0'  STOP 

GO  TO  25 

END 


FUNCTION  F ' X j 
COMMON  A.B 

V  -  . 5*  ( A+8+ '  B-A  >*X  ) 

C  ENTER  INTEGRAND  NEXT  AS  FCN.  OF  V 
G  •  2./C2.+SINC  10.*AC0S<-1.  UVD 
F  •  GISQRT ' 1 . -XT*2  ) 

RETURN 

END 


R;  T*0.02/0.19  14M3:34 


Exhibit  1 


OKI  PA 5*  IS  8ASI  iJUALil  1  : 
TBdli  SttPI  fWUkiSiiSD  TV  DDC 
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TION  subprogram.  This  feature  may  cause  the  sequence  of  Gauss- 
Chebyshev  quadratures  to  converge  more  slowly  than  the  correspond¬ 
ing  sequence  of  Gauss-Legendre  quadratures  discussed  in  the  next 
section  (where  the  weight  function  =  1) ,  but  the  simplicity  of 
formula  (3.2)  still  is  a  cogent  argument  for  its  use. 

Table  1  gives  some  sample  output  for  an  assortment  of  defin¬ 
ite  integrals.  Since  the  computation  was  carried  out  in  single 
precision  IBM  370  arithmetic,  only  5  significant  figures  are 
given.  Observe  that  in  all  cases  except  the  last  convergence 
(to  5  figures)  is  achieved  by  an  order  200  formula.  The  time 
column  gives  the  execution  times  required  to  compute  the  last 
entry  in  each  row.  These  times  are  a  bit  higher  than  necessary 
in  most  cases  because  of  the  large  increases  in  N.  That  is,  if 
we  increased  N  more  slowly,  say  by  increments  of  5  or  10,  we 
would  generally  achieve  convergence  to  5  places  sooner  than  shown 
in  Table  1,  with  a  correspondingly  reduced  execution  time. 


EXAMPLES  OF  GAUSS-CHEBYSHEV  QUADRATURE 


(l-jyi)dy  .99161  1.00263  .99967  .99991  .99995  .99997  1.0000  .012  sec 


w 


4.2  Unilateral  Error  Bounds 

A  final  remark  about  Gauss-Chebyshev  quadrature  pertains  to 
the  possibility  of  a  one-sided  error  bound,  in  the  sense  that  for 
some  integrands  f  it  may  happen  that  either 

G  (f;  Chb. )  <  f  f  fo)(?X  ,  (3.3) 

n  '/-lsjl-x2 

or 

G  (f;  Chb.)  >  f1  . 

n  ~  j-l  ,  (3.4) 

for  all  n.  For  example,  inequality  (3.4)  is  valid  for  the  first 
three  integrals  in  Table  1,  and  also  holds  for  the  fourth  integral, 
provided  that  n>4, 

These  inequalities  would  obviously  follow  from  the  stronger 
results  of  monotone  convergence 

G  (f;  Chb.)  <  G  ..  (f ,*  Chb.)  (3.3m) 

n  —  n+l 


or 


Gn(f?  Chb.)  >Gn+1(f;  Chb.),  (3.4m) 

for  all  n  (or,  for  all  sufficiently  large  n) .  In  this  section  we 
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give  some  conditions  of  f  sufficient  to  ensure  one  of  the  pre¬ 
ceding  inequalities  (3.3)  or  (3.4). 

In  general,  let 

EG  (f;  Chb. )  =  f  f  (x)  dX  -  G  (f;  Chb.) 
n  -1  1-x2  n 

be  the  error  in  the  Gauss-Chebyshev  rule.  One  expression  for 
EGn(f;  Chb.)  was  given  in  (3.2),  as  a  special  case  of  (2.6).  For 
the  general  case  of  Gauss-Jacobi  quadrature  (where  the  weight 
function  <o  (x)  =  (l-x)a  (1+x)^,  -l<a , B , | x | £l ) ,  we  can  deduce  (2.6) 
and  hence  (3.2)  from  the  Peano  kernel  theorem  [1,  p.  18]  which 
permits  us  to  write 


EG  (f ;  cj  ) 
n 


f  (2n) 


( t)  dt 


(3.5) 


for  every  fe  C2n[-1,  1],  Here  is  a  certain  "kernel"  which 

depends  on  the  weight  fucntion  to  and  the  order  n.  It  is  known 
that  is  of  class  C2n_2  on  [-1,  1],  is  positive  for  [  1 1  <  1 ,  and 

U) 

in  the  a=B  (ultraspherical)  case,  '  is  an  even  function  [12, 

Ch.  4]. 

Now  consider  the  special  integrand  f(x)  =  x^,  p=0,  1,  2,***. 
Applying  (3.5)  we  have 


EG  (xP;  to)  =  p(p-l)  *•*  (p-2n+l)  f  K  (n)  (t)  p-2ndt 


n 


-1 
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2p (p-1) • • • (p-2n+l)  /  K' ;  '(t)tF  dt  >0,  if  p  is  even, 

•'o 

0  ,  if  p  is  odd, 

provided  that  is  an  even  function.  These  equations  remain 

true  if  f  is  replaced  by  any  polynomial  with  positive  coefficients 
or  a  limit  of  such  polynomials.  Hence,  for  such  functions,  in¬ 
equality  (3.3)  is  valid  (note  that  this  result  has  been  demon¬ 
strated  to  hold  for  all  Gauss  quadrature  formulas  derived  from 
ultraspherical  polynomials,  in  particular  for  the  Gauss-Legendre 
rule  of  the  next  section) . 

A  general  class  of  integrands  for  which  (3.3)  is  valid  is, 
for  instance,  the  class  of  analytic  functions  f  on  (-1,1)  with 

f (p) (0)  ^  0,  p  =  1,  2,  3,...  .  (3.6) 

And  by  starting  with  f(x)  =  -xP,  p  =  1,  2,  3,...,  and  repeating 
the  above  analysis,  we  can  obtain  a  class  of  functions  for  which 
inequality  (3.4)  is  valid,  namely,  the  class  of  analytic  functions 
f  for  which  (3.6)  holds  with  the  opposite  inequality  sign. 
Empirical  evidence  suggests  that  in  fact  the  stronger  inequalities 
(3.3m)  and  (3.4m)  hold  for  these  two  classes  of  integrands  res¬ 
pectively. 
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Another  condition  sufficient  to  ensure  one  of  the  estimates 

(3.3)  or  (3.4)  follows  immediately  from  the  error  formula  (3.2) . 

( 2n ) 

Namely,  if  f  has  derivatives  of  all  orders  on  (-1,  1)  and  f  (x) 

> _  0  (resp.  £  0)  on  (-1,  1),  then  inequality  (3.3)  (resp.  (3.4)) 
is  valid.  Indeed  this  condition  may  be  established  for  quite 
general  weight  functions,  as  was  originally  shown  by  Shohat;  see 
[5,  p.  92]. 
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V 


GAUSS-LEGENDRE  QUADRATURE 


We  next  consider  integration  with  respect  to  the  simplest 
weight  function  on  [-1,  1],  namely  w(x)  =  1.  In  this  case  the 
resulting  quadrature  formula  (2.5)  will  approximate  the  integral 


f (x) dx 


(4.1) 


The  corresponding  orthonormal  polynomials  associated  with  this 
weight  function  are  the  Legendre  polynomials.  They  can  be  ex¬ 
pressed  as 


P  (x) 
*n 


p  (x)  , 
n 


(4.2) 


where 


Pn(x) 


dx11 


(x2-l) 


n 

r 


(4.3) 


this  last  expression  being  the  Rodrigues  formula  [6,  v.l,  Ch.  15]. 

In  order  to  effectively  compute  the  abscissas  and  weights  of 
the  corresponding  quadrature  formula  (2.5) ,  we  must  utilize  the 
eigenanalysis  of  part  II.  This  means  that  we  need  the  recurrence 
relation  (2.2)  from  which  we  can  construct  the  Jacobi  matrix  A  in 
(2.4).  Now  the  recurrence  for  the  polynomials  Pn  in  (4.3)  is 


nP  (x)  =  (2n-l)  xP  , (x)  -  (n-l)P  ,  (x) 
n  n-l  n— 2 


[6],  Combining  this  with  the  expression  for  the  normalized 
polynomials  pn  in  (4.2)  yields  the  desired  recurrence  relation: 


Pn(x) 


2n+l 

2n-l 


xpn-i<xl 


2n+l 

2n-3 


Pn.2(X>- 


(4.4) 


From  (4.4)  follows  the  consistency  relation  (2.3)  and  from  this 
follows 


_n  =  1  =  n-l 

an  an-l  ^(2n-l)  (2n-3)  * 


Since  the  {bn>  of  the  general  recurrence  relation  (2.2)  are  zero 
here,  we  have  all  the  entries  available  to  construct  the  desired 
matrix  A. 


5.1  Program  and  Examples 

Exhibit  2  gives  a  FORTRAN  program  which  enables  a  user  to 
calculate  an  n*"*1  order  Gauss-Legendre  approximation  to  an  integral 
of  the  form 
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FUNCTION  FIX> 

IMPLICIT  «E»Lt8<A-H.O-Z> 
CONNCW  A,  I 


There  are  two  subroutines:  subroutine  FORM  calculates  the  weights 
(stored  in  the  array  WGTS)  and  the  abscissas  (stored  in  DD) ; 
the  IMSL  subroutine  EQRT2S  is  called  by  FORM  and  performs  the 

required  eigenanalysis  on  the  matrix  A.  The  program  is  written 

in  double  precision  because  of  the  added  complexity  of  the  eigen- 
routine. 

Table  2  gives  the  results  of  using  this  program  to  compute 

the  same  integrals  as  given  in  Table  1.  The  time  column  gives 

the  execution  time  required  for  the  last  value  reported  in  each 
row. 


5.2  Comparison  with  Gauss-Chebyshev  Quadrature 

Naturally,  Tables  1  and  2  invite  comparison.  Although  the 
evidence  is  limited  we  draw  the  following  general  conclusions: 

a)  Because  of  the  artifical  weight  function  introduced  in 
Gauss-Chebyshev  formula,  its  convergence  as  measured  by 
the  order  of  the  formula  (and  hence  the  number  of  func¬ 
tion  evaluations)  is  much  slower  than  that  of  the  Gauss- 
Legendre  formula.  The  latter  however  requires  added 
calculation  to  establish  the  weights  and  abscissas. 

Thus  we  need  to  decide  (in  advance!)  whether  this  add¬ 
itional  calculation  is  worthwhile.  This  depends  on  the 
integrand.  As  the  complexity  of  its  evaluation  at 
various  points  increases,  so  does  the  gain  in  using 
fewer  function  evaluations,  hence  the  Gauss-Legendre 
formula. 
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b)  Also  if  the  integrand  is  analytic  on  a  neighborhood  of 
the  segment  [a,  b]  (as  are  the  second  and  third  in¬ 
tegrands  in  the  Tables) ,  then  the  rapid  convergence  of 
Gauss-Legendre  argues  in  its  favor.  These  two  ex¬ 
amples  show  time  ratios  of  about  3:1  in  favor  of  Gauss- 
Legendre  . 

c)  When  the  integrand  is  analytic  on  (a,  b)  but  not  on 
[a,  b]  (as  are  the  first  and  fourth  integrands  in  the 
Tables) ,  the  advantage  is  still  with  Gauss-Legendre 
but  somewhat  more  narrowly. 

d)  When  the  integrand  is  badly  behaved  with  respect  to 
polynomial  approximation  (as  are  the  final  two  in¬ 
tegrands)  ,  the  Gauss-Chebyshev  method  becomes  much 
preferred,  because  of  the  ease  with  which  its  order 
is  increased. 

The  general  error  formula  (2.6)  for  Gauss-Legendre  quadrature 


dx 


Gn(f;  Leg.) 


+  2*4n» (n!)4  f (2n) (g) 
( (2n) ! ) 3  (2n+l) 


(4.5) 


for  1 5 1  <1  and  feC2n[-l,  1],  and  this  too  invites  comparison  with 
the  corresponding  formula  (3.2)  for  Gauss-Chebyshev  quadrature. 

Let  GCR  (resp.,  GLR)  be  the  coefficient  of  f^2n^  (£)  in  (3.2) 

(resp. ,  in  (4.5)),  and  consider  the  ratio  GCR/GLR.  Making  use 
of  Stirling's  approximation  it  is  easy  to  see  that 


i-t 
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The  availability  of  subroutine  FORM  in  the  program  of  Exhibit 
2  makes  it  very  easy  to  implement  the  product  Gauss  formula  of 
(2.11)  for  the  case  oj  (x)  =  m  (y)  =  1,  and  the  rectangle  a<y<b, 
c<x<d.  Similarly  one  could  implement  a  product  Gauss-Chebyshev 
rule,  although  we  have  not  done  so  for  this  report,  partly  from 
a  reluctance  to  artificially  introduce  two  weight  functions,  and 
partly  because  of  the  increased  number  of  function  evaluations 
involved  in  high  order  product  integral  formulas.  But  as  before, 
if  too  high  an  order  of  Gauss-Legendre  is  required  for  the  desired 
accuracy,  then  the  advantage  in  computation  time  will  pass  to 
Gauss-Chebyshev . 

There  are  several  alternate  methods  of  constructing  quadrature 
formulas  for  multiple  integrals  [10].  Three  popular  ones  are 
a)  choosing  the  weights  and  nodes  (=  points  at  which  the  integrand 
is  evaluated)  so  as  to  make  the  resulting  formula  exactly  integrate 
polynomials  of  a  certain  degree  or  less;  b)  minimizing  the  max- 
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imum  error  over  a  certain  class  of  possible  integrands;  and  c) 
statistical  (Monte  Carlo)  sampling. 

We  looked  at  several  methods  of  type  a) ,  as  their  theory 
is  the  simplest  (of  the  three  general  methods) ,  and  closest  in 
spirit  to  Gaussian  quadrature  with  its  greatest  success  for  in¬ 
tegrands  having  good  polynomial  approximation  properties.  Pro¬ 
cedures  for  constructing  quadrature  rules  of  the  form 


r 

J  f  (u)  du  =  f  (P  ) 

r*  '  _ T  -1*  1 


(4.6) 


which  are  exact  for  all  polynomials  of  degree  _<  d,  and  for  which 
N  is  a  minimum  are  quite  involved,  and  in  fact  have  only  been 
carried  out  in  a  few  special  cases  [10].  Particularly  interesting 
from  the  viewpoints  of  accuracy  and  simplicity  are  formulas  due 
to  Radon  [10,  §  3.12]  and  Albrecht-Collatz  [10,  §  8.2],  These 
are  7  point  (N=7  in  (4.6))  formulas  of  degree  5  for  the  unit  square 
-l<x,  y£l.  That  is,  they  will  integrate  any  quintic  polynomial 
in  (x,  y)  exactly  over  this  square,  and  further,  no  formula  of 
type  (4.6)  with  N<7  can  do  this. 

The  Albrecht-Collatz  formula  is 


//fU,  y )  dxdy  a  4  |f(0,  0)  +  ^(f(r,  r)  + 


-1  -1 


f(-r,  -r))  +  jg-(f(s,  -t)  +  f(-s,  t)  +  f(t,  -s)  +  f(-t,  s)  )  , 
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where  r2  =  7/15,  s2  =  (7+\kI)/15,  t2  =  (7-\^T)/15. 

In  Table  3  we  present  results  from  using  the  product  Gauss- 
Legendre  and  the  Albrecht-Collatz  formulas  on  a  few  double  in¬ 
tegrals.  As  expected,  the  convergence  is  rapid  for  analytic  in¬ 
tegrands,  and  the  Albrecht-Collatz  formula  usually  gives  an 
accuracy  of  +  .05%  for  such  integrands,  although  there  is  no 
practical  way  to  predict  its  accuracy  ex  ante.  The  fourth  integral 
exhibits  quite  slow  convergence  due  to  the  pole  at  (1,  1) .  The 
fifth  integral  is  typical  of  the  sort  appearing  in  optical  de¬ 
tection  theory.  It  represents  the  output  of  a  rectangular  de¬ 
tector;  the  integrand  is  the  irradiance  function,  is  the  Bessel 
function  of  first  order,  and  e  is  the  obscuration  factor  [4], 


EXAMPLES  OF  TWO-DIMENSIONAL  GAUSS-LEGENDRE  QUADRATURE 


dxdy  .44314  ,47133  .48779  .81605  .52697  .53155  .33136  .53333 


VI .  GAUSS-LAGUERRE  QUADRATURE 


In  this  final  section  we  consider  integration  over  the  half¬ 
line  [a,  °°)  with  the  respect  to  the  weight  function  w(x;  a)  = 
a  -x 

x  e  .  The  constraint  ot>-l  is  enforced  as  otherwise  co  ( • ;  a)  is  not 
integrable  over  [a,  °°)  . 

6.1  Special  Features 

There  are  two  new  features  in  the  present  situation:  the 
infinite  extent  of  the  region  of  integration,  and  the  parameter 
a.  As  already  observed  in  paragraph  2.4,  the  first  feature  creates 
difficulties  in  establishing  the  convergence 

00 

lim  G  (f;  oj  ( • ;  a))  =  J  f  (x)w(x;  a)dx, 
n-*»  n  a 


but  these  difficulties  are  of  a  theoretical  nature  and,  in  any 
event,  have  been  resolved,  as  noted  earlier. 

The  use  of  the  free  parameter  a  is  of  greater  practical  in¬ 
terest.  We  have  not  explicitly  encountered  its  analog  before,  al¬ 
though  one  is  present.  Namely,  the  Chebyshev  and  Lengendre  poly¬ 
nomials,  already  used,  are  orthogonal  over  [-1,1]  with  respect  to 

the  weight  functions  (1-x  )  and  1,  respectively.  Hence  they  are 
special  cases  of  polynomials  orthogonal  over  [-1,  1]  with  respect 


2  ^ 

to  the  weight  function  (1-x  )  ,  B>-1.  Any  such  set  of  polynomials 
is  called  a  set  of  ultraspherical  polynomials,  and  in  particular 


F 
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obeys  a  recurrence  relation  of  the  form  (2.2),  from  which  a  corres¬ 
ponding  Gaussian  quadrature  formula  could  be  developed  as  has  al¬ 
ready  been  done  for  the  special  cases  Q=-h  and  g=0.  The  only 
practical  effect  of  doing  so  is  to  slightly  shift  the  corresponding 
abscissas  inside  the  interval  [-1,  1],  (See  [9,  p.  121]  for  the 
pertinent  inequalities  relating  these  abscissas  to  those  of  the 
Chebyshev  polynomials  of  the  first  and  second  kinds.  These 
abscissas  are  always  symmetric  (with  respect  to  0)  for  the  ultra- 
spherical  polynomials.)  Unless  the  weight  function  (l-x^)  is 
explicitly  present  in  an  integrand,  there  seems  little  point  to 
introducing  it  artificially  (as  emphasized  earlier,  the  case  B=-Sj 
presents  a  striking  exception  to  this  advice) . 

A  discussion  of  the  zeros  of  the  n1"  generalized  Laguerre 

polynomial  occurs  in  [9,  §  6.31],  In  fact,  if  x/n^  is  the 

n  k 

ktk  zero  of  then  the  following  bounds  are  known: 

(jv/2)2  {n) 

n  +  (a  +  l)/2  <  xk  K 


(k+(a+l)/2)  •  ■2k+at1.'>'. (.(2k+a+l)  +h-a  )_  , 

n+ (a+1) /2 

where  {j^.}  is  the  sequence  of  positive  zeros  of  the  Bessel  function 
J  .  As  is  also  known  [9,  §  6.21],  is  an  increasing  function 

CX  K 

of  a. 
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( n ) 

This  last  property  of  x^  ;  is  of  practical  importance  for 
Gauss-Laguerre  quadrature.  For  a  fixed  order  n  the  location  of 


the  abscissas  can  be  shifted  towards  (or  away  from)  the  left- 
hand  end  point  by  letting  a  decrease  (resp.,  increase).  In  par¬ 
ticular  the  location  of  the  abscissas  can  be  chosen  to  reasonably 
coincide  with  the  range  of  greatest  variation  of  the  integrand. 

To  illustrate  this  behavior  of  abscissas  (and  weights)  of 
Gauss-Laguerre  quadrature  rules,  we  display  in  Table  4  the 
abscissas  and  weights  for  several  different  8  point  formulas. 

This  table  shows  how  the  abscissas  drift  off  to  the  right  as  a 
increases,  and  also  how  the  corresponding  integrand  values  are 
given  higher  weight.  The  right  hand  column  gives  the  theoretical 
sum  of  the  weights,  computed  from  the  formula 

n  jo 

S  =  f  xae  xdx  =  T  (a+1)  . 

k=l  k  */0 

In  all  cases  the  actual  sum  of  weights  agrees  with  the  theoretical 
sum  to  within  1  or  2  units  in  the  6*"^  decimal  place.  We  can  see 
the  sum  of  the  weights  first  decreasing  to  a  minimum  (actual  value 
=  0.885603  when  a  =  0.461632),  and  then  increasing  without  bound. 
This  information  follows  from  standard  properties  of  the  gamma 
function  [13,  Ch.  6], 
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6.2  Program  and  Examples 

Table  4  was  derived  from  the  FORTRAN  program  in  Exhibit  3. 
This  program  enables  a  user  to  compute  an  approximation  to  the 
integral 


the  user  must  supply  as  input  a  value  of  the  parameter  a ,  the 
order  N  of  the  corresponding  Gauss-Laguerre  formula,  and  the 
value  T  of  the  lower  endpoint.  As  in  the  program  of  Exhibit  2 
for  Gauss-Legendre  quadrature ,  there  are  two  subroutines :  FORM 
and  EQRT2S  (from  IMSL) ,  which  together  effect  the  computation 
of  the  relevant  weights  and  abscissas. 

Table  5  displays  a  few  examples  of  Gauss-Laguerre  quadrature. 
The  first  integrand  is  a  nice  entire  function  for  which  the  choice 
a  =  0  leads  to  rapid  convergence.  The  second  integrand  exhibits 

slowly  decaying  oscillatory  behavior  and  convergence  is  poor. 
Probably  a  =  0  would  be  the  default  choice  here  too,  although  the 
integrand  is  inherently  badly  behaved  for  quadrature.  The  next 
integrand  has  a  singularity  at  x  =  0;  the  feature  in  general  would 
suggest  a  value  of  a  <  0,  and  the  factor  x  in  the  integrand  in 
turn  leads  to  a  =  -.5  as  the  preferred  choice.  The  last  integrand 
is  (up  to  a  constant  factor)  the  type  occuring  in  integration  of 


<r 
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EXAMPLES  OF  GAUSS-LAGUERRE  QUADRATURE 


Planck's  equation  for  radiant  graybody  sterance.  The  interval  of 

integration  [2,°°)  corresponds  roughly  to  the  wavelength  band 

(0,  7194].  Surprisingly  enough  here,  the  "natural"  choice  a  =  3 

is  considerably  less  effective  than  smaller  a's.  The  reason  for 

3  x 

this  can  be  deduced  from  Table  4;  since  the  integrand  x  /(e  -1) 
is  rapidly  decreasing  after  its  peak  at  x  *  2.8,  it  is  advantageous 
to  have  the  abscissas  clustered  closer  to  the  left  end  point 
(=2,  in  our  example) ,  and  as  we  know,  this  requires  smaller 
values  of  a. 
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VII.  FINAL  EXAMPLES  AND  SUMMARY 

In  this  report  we  have  advanced  a  fast  and  efficient  method 
for  carrying  out  numerical  integration  by  means  of  Gaussian 
quadrature.  We  have  stressed  that  Gaussian  methods  should  only 
be  applied  in  situation  where  the  integrand  has  the  property  of 
good  polynomial  approximation.  In  the  contrary  cases,  when  either 
the  integrand  is  insufficiently  differentiable  or  has  a  singular¬ 
ity  on  or  near  the  boundary  of  the  region  of  integation,  other 
methods  should  be  adopted.  In  particular,  in  the  first  case  some 
version  of  the  trapazoid  rule,  such  as  Romberg  integration,  might 
be  tried  (see  Introduction) .  On  the  other  hand,  the  treatment  of 
singularities  is  often  a  matter  of  ad  hoc  techniques,  several  of 
which  are  suggested  in  [14], 

Sometimes  a  simple  change  of  variable  can  greatly  expedite 
the  quadrature  process.  We  give  next  two  final  examples,  both 
of  which  can  be  viewed  either  as  integrands  with  singularity  at 
one  end  point  of  a  finite  interval  or  as  integrands  analytic  over 
an  infinite  interval.  In  the  first  case  the  change  of  variables 
does  not  help  at  all,  while  in  the  second  case  the  improvement  is 
dramatic. 
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Example  1 

/•l  i  -co  i 

J  (l-e"x-e"  x)^£  =  f  (l-e-X-e“  x)^  =  Y  =  .577216, 

*'o  x  J1  x 

where  y  is  Euler's  constant.  The  results  of  three  attempts  at 
integration  by  our  methods  are  given  next. 

Gauss-Chebyshev 


N= 

10 

25 

50 

75 

100 

150  | 

1  200 

Value 

.579836 

.577631 

.577316 

•577257 

i 

.577235 

.577221 

.577214 

(time=,024  sec.) 

Gauss-Legendre 


N= 

5 

6 

7 

8 

9 

10 

11 

12  i 

15 

Value 

.577684 

.577450 

.577155 

.577165 

.577218 

( 

.5772277 

.577218 

i 

- f 

.577214- 

.577216 

(time=.055  ) 
sec. 
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Gauss-Laguerre 


N= 


I  I 

I 

8  i  10  12 


15 


20 


I 


Value 


a 


-.9 

-.75 

-.5 

-.25 

0. 


.52544 

! 

.52999 

1  .53823 

.53134 

|  .53923 

.54485 

.52840 

.52257 

1 

• 

.53168 

i 

.55260 

I 

e54848  1 

.55456 

t 

*55855 

.54907 

1  .55494 

.55470 

.  J  J  jZu 

'  .55009 

1 

.56126 

.56249 

.56267 

.56149 

.55931 


.56545 

.56635 

.56646 

.56555 

.56389 


We  can  see  from  these  results  that  conversion  from  the  in¬ 
terval  (0,1]  to  [1,°°)  has  not  achieved  any  worthwhile  results; 
the  Gauss-Laguerre  formulas  are  all  converging  very  slowly.  The 
Gauss-Chebyshev  method  has  essentially  converged  by  N  =  200;  the 
error  in  the  6  decimal  place  is  due  to  the  inaccuracy  of  the 
single  precision  arithmetic. 


Example  2 

f1  ( log  x)ndx  =  (-1 )n  f  yne"ydy  =  (-l)nnl 

*/0 

In  this  example  we  see  that  the  change  of  variable  converts  an 
intractable  integral  over  (0,1)  into  a  very  simple  integral  over 
[O,00)  which  can  either  be  evaluated  in  closed  form  (using  integra- 


In  conclusion,  if  asked  to  recommend  a  single  all-purpose  integra¬ 
tion  routine,  our  choice  would  have  to  be  the  Gauss— Chebyshev  for¬ 
mula  (even  for  infinite  intervals,  by  truncation).  Its  simplicity 
and  the  speed  with  which  its  order  can  be  increased  are  the  de¬ 
cisive  factors. 
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